Modeling Extreme Mass Ratio Inspirals within the Effective- One- Body Approach 
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We present the first models of extreme-mass-ratio inspirals within the effective-one-body (EOB) 
formalism, focusing on quasi-circular orbits into non-rotating black holes. We show that the phase 
difference and (Newtonian normalized) amplitude difference between analytical EOB and numerical 
Teukolsky-based gravitational waveforms can be reduced to < 10 _1 rads and < 2x 10 -3 , respectively, 
after a 2-year evolution. The inclusion of post-Newtonian self-force terms in the EOB approach leads 
to a phase disagreement of ~ 6-27 rads after a 2-year evolution. Such inclusion could also allow for 
the EOB modeling of waveforms from intermediate-mass ratio, quasi-circular inspirals. 



Introduction. Extreme mass ratio inspirals (EMRIs) 
are the gravitational wave (GW) driven coalescences of 
stellar mass compact objects with supermassive black 
holes (SMBHs). When the large black hole's (BH's) 
mass is in the range 10 5 M© — 10 7 Mq, EMPJ waves are 
emitted at frequencies well suited to measurement by the 
planned Laser Interferometer Space Antenna (LISA). Be- 
cause EMRI events are expected to be abundant [l| and 
will carry detailed information about strong-field space- 
times near SMBHs 0, they are high-priority targets for 
LISA observation. The intrinsic feebleness of these waves 
will require accurate waveform templates to detect and 
faithfully measure the signals produced by Nature. Be- 
cause EMRIs can spend thousands of cycles in the close 
vicinity of the SMBH's innermost stable circular orbit 
(ISCO), traditional post-Newtonian expansions are not 
well suited to modeling their signals; the binary's orbital 
speed is v/c ~ 0.1-0.5, a regime where traditional post- 
Newtonian techniques perform poorly. Numerical models 
built using BH perturbation theory should be able to re- 
liably model EMRI signals. However, the computational 
cost of covering the full span of EMRI parameter space 
(including effects due to BH spin, non-equatorial orbits, 
and eccentricity) is likely to be very high [![. 

This has motivated us to examine techniques for reli- 
ably approximating these waves at a much smaller com- 
putational cost. The effective-one-body (EOB) formal- 
ism was introduced as a way to analytically describe 
the inspiral, merger, and ringdown waves emitted by 
comparable-mass BH binaries 0, This formalism 
was then extended to higher post-Newtonian (PN) or- 
ders [5J1, spinning BHs [a-lSj, small mass-ratio merg- 
ers P.Tlol| , and was further improved by introducing fac- 
torized waveforms 0, 11 1. By calibrating a few ad- 
justable parameters in the EOB-dynamics and wave- 
forms, [12, EH showed that the phase and amplitude 
of the EOB and numerical-relativity waveforms can be 
made to agree within the numerical error of the simu- 



lations, thus providing GW detectors with faithful tem- 
plates. In this analysis, we consider calibrating EOB 
with BH perturbation theory templates in order to simi- 
larly model EMRI waves. Such an analysis must be done 
separately from the previous EOB-numerical relativity 
calibration, because sufficiently long numerical relativity 
simulations are currently not available for EMRIs. This 
is because the large mass ratio in EMRI systems leads to 
tens to hundreds of thousands of detectable cycles within 
a one- year window using LISA, thus requiring extremely 
long simulations that are currently computationally pro- 
hibitive. 

As a first step, we restrict our models to a small com- 
pact object spiraling along a quasi-circular orbit into a 
non-spinning SMBH 14|. Although the assumptions of 
circularity and zero spin can and will be relaxed in the 
future, there exist astrophysical motivations for this ini- 
tial choice of binary configuration. For example, the tidal 
separation scenario for EMRIs [lj| implies nearly circular 
but arbitrarily inclined orbits in the > 10~ 4 Hz frequency 
band relevant for LISA, and the accretion disk capture 
picture implies orbits that are both nearly circu- 

lar and in the equatorial plane of the SMBH. In addition, 
the characteristics of the SMBHs themselves are uncer- 
tain in the < 10 7 Mq mass range most relevant for LISA. 
In some astrophysical scenarios, the growth of these BHs 
is dominated by the accretion of stars moving on ran- 
dom trajectories, instead of by the accretion of gas disks, 
thought to be more important for higher-mass SMBHs 
[Hj]. Such growth would lead to a = \J\/M 2 < 1 (in 
natural units with G = c = 1, which we use through- 
out this paper), hence the Schwarzschild (nonrotating) 
spacetime is a reasonable first approximation. 

We now systematically compare EMRI waveforms 
computed in the EOB approach to those calculated us- 
ing BH perturbation theory via the numerical solution of 
the Teukolsky equation [20h23I|. As we describe in the 
remainder of this Letter, we find that appropriately cali- 
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brated EOB waveforms in fact do an excellent job model- 
ing waves computed using BH perturbation theory. This 
suggests that the EOB scheme is very likely to be an out- 
standing tool for modeling EMRI waves in future LISA 
data analysis. 

Analytical and Numerical Modeling. For a BH bi- 
nary with masses mi and 777.2, we set M = mi + m% 
and fj, — m\m%jM. — vM. In absence of spins, 
the motion is constrained to a plane. Let us intro- 
duce Schwarzschild-like coordinates (r, (where r is 
M -normalized) centered on the binary's center of mass, 
as well as their reduced (/x-normalized) conjugate mo- 
menta (p r ,P§)- The non-spinning EOB Hamiltonian 



then reads |3j H Ical = M y/l + 2v\(H eB - Jtj/Ji] 
where the effective Hamiltonian is [H, H, 



A(r) 



1 + ^ + 2(4-3,),^ 



M, 



(1) 



We use here the reduced conjugate momentum p r ^ to the 
EOB tortoise radial coordinate r* because it improves the 
numerical stability of the code [l(|. The tortoise coor- 
dinate is defined via dr„/dr = y/D(r)/A(r), where A(r) 
and D(r) are obtained by applying the Pade resumma- 
tion [5] to the Taylor-expanded forms 0, [B| 



A T (r) = 1 - - + — 
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D T (r) = 1-^ + 2^(3^-26)4 



(2) 
(3) 



The EOB Hamilton equations are written in terms of 
the reduced (dimensionless) quantities jj rcal = H rcld /fi, 
t = t/M 0: 



dr 
dt 

dPr, 

dt 



real 



A{r) dH 
v/^W dp r , ' 

A(r) dH Ical 
~ ^D{r) dr 



dt 



dH 



real 



dp$ 
dt 



dp® 



(4) 
(5) 



where J-$ is a Pade-resummed radiation-reaction force 0, 
24j . related to the GW energy dissipation to be defined 
later. Initial data is obtained through a mock evolu- 
tion, which is initialized at an initial orbital separation 
of 100M using initial conditions for a quasi-circular in- 
spiral 

With the EOB inspiral dynamics in hand, we com- 
pute the multipole-decomposed GW hg m (I and m re- 
fer to spherical harmonics), following the factorized PN 
prescription of [ll|, which depends directly on orbital 
quantities. The EOB GW phase is computed by solving 
&£m = — (l/ m ) ^AP-tm/him]- Errors in the EOB wave- 
forms arise due to inaccuracies in the numerical solution 
of Eqs. (jH) and (0 and inaccurate initial data. We have 
investigated such sources of error and estimate them to 
be no worse than 5<&22 % 0.03 rads in the waveform's 



phase and <5/i22/ft-22 ^ 10~ 7 in the normalized amplitude 
after a 2-year evolution. This cumulative error is pri- 
marily dominated by the accuracy of the routine used in 
Mathematica to solve Eqs. (|4]) and (0. 

We compare EOB waves with waveforms computed 
in BH perturbation theory by solving the Teukolsky 
equation. We use the code described in [22} (modified 
with the spectral techniques of 0) to construct the 
Newman- Penrose curvature scalar ^4. Our code works 
in the frequency domain, decomposing ^4 into ^4 — 



R 1 J2i m z - 



iral 



where -?X 



a spin-weight —2 spherical harmonic, f2 is the frequency 
of circular Schwarzschild orbits and R is the distance 
from the center of mass to the observer. The ampli- 
tude Zi m is found by first constructing a Green's func- 
tion to the radial Teukolsky equation, and then integrat- 
ing that Green's function over a source made from the 
stress-energy tensor of the small body orbiting the BH; 
see 



22] for specifics. 



The radial Teukolsky equation possesses two asymp- 
totic solutions that determine the behavior of "04 at spa- 
tial infinity and near the event horizon. In the distant 
radiation zone, ^4 is related to the GWs carried away 
from the system via -04 — > l/2(/i + — ih x ). Therefore, the 
solution to the radial Teukolsky equation that describes 
purely outgoing radiation at spatial infinity can be used 
to construct the flux of radiation and the waveform that 
distant observers measure. On the other hand, near the 
event horizon, ^4 describes tidal interactions of the BH 
with the orbiting body 2l| . Thus, the solution to the 



radial Teukolsky equation that describes purely ingoing 
radiation at the horizon can be used to construct the ra- 
diation flux absorbed by the BH. With these fluxes, we 
can then calculate the rate at which the orbital radius 
changes, r, by noting that for slow backreaction the sys- 
tem evolves through a sequence of geodesic orbits. 

We construct the ip4 solution on a discrete grid of orbits 
from r = 10,000M to the Schwarzschild ISCO at r = 
6M (in Boyer-Lindquist coordinates), evenly spacing our 
orbits in v = y/M/r. (Since stable circular orbits do not 
exist for r < 6M, we cannot infer r from dE/dt in this 
regime.) Errors in the Teukolsky-based waveforms are 
dominated by truncation of the (£, m) sums in ^4 and 
due to the discretization of the orbital phase space when 
the fluxes are cubic-spline interpolated from a discrete 
adiabatic sequence of geodesic orbits. The sums and the 
discretization are chosen such that the fractional error in 
the flux is smaller than 10 -10 22, [23|. In practice, in 
the low velocity region v < 0.1, we find that the flux is 
accurate to at least 10 -13 . Such an error translates to 
inaccuracies in the GW phase of less than 10~ 2 rads over 
a 2-year evolution. 

Systems, Regions and Models. To demonstrate the 
flexibility of the EOB model in matching the Teukolsky- 
based waveforms, we examine two fiducial EMRI sys- 
tems, labeled system-I and system-II, that sample dif- 
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1: Absolute value of the difference in the Newtonian 
normalized Teukolsky and EOB fluxes as a function of or- 
bital velocity. Calibrating the Pade or p— flux improves the 
agreement by orders of magnitude 

ferent regions of the LISA noise curve. In both cases, 
we consider a 2-year long quasi-circular inspiral of non- 
spinning BHs. System-I has (mi, ma) = (10 5 , 10)M Q ; 
system-II has (mi, m^) — (10 6 ,10)M Q . We do not con- 
sider lower or higher total mass binaries as they would ci- 
ther reach the ISCO outside the LISA optimal sensitivity 
band (LISA's noise rises sharply above ~ 10 -2 Hz) or lie 
significantly inside the white-dwarf confusion limit (much 
below ~ 0.002 Hz Q). System-II {m 2 /m 1 = 10~ 5 ) 
begins at an initial separation r; n ~ 10. 6M and ter- 
minates at the ISCO, sweeping GW frequencies in the 
range / Gw € [1.8 x 10" 3 ,4.4 x 10" 3 ] Hz. System-I 
(m2/mi = 10~ 4 ) starts at r- m ~ 29.34A/ and termi- 
nates at rfi n ~ 16. 1M, sweeping frequencies in the range 
/gw £ [4 x 10~ 3 , 10~ 2 ] Hz. The mass ratios we consider, 
(10~ 4 ,10~ 5 ) are two orders of magnitude smaller than 
those studied in the complementary analyses of 0, [l(| ■ 
As such, our in-band signal is dominated by a long in- 
spiral; the contributions of the final plunge, merger, and 
ringdown, which dominate the signal of [9|, [l0( , are much 
less important here. These choices allow the study of the 
early and late EMRI dynamics, while guaranteeing the 
GW signal is in the sensitive part of the LISA band. 

We define two EOB models differing in the resumma- 
tion of the radiation-reaction force in Eq. ([5]). Using 
the balance law, we write J-$ = —F/(M), where F is 
the GW energy flux. We use (i) the Pade-approximant 
to the energy flux [13j, |24j F p = FP(v po i c ), where w po i 
is an adjustable parameter locating the EOB light-ring, 
and p + q is twice the approximant's PN order [i.e., 
(v /c) (P +q > ] , and (ii) the p-approximant to the energy 
flux 0P = 2/(16tt) Y?^lY2zi(™n) 2 \Rh em \ 2 . Ex- 
cept when investigating the effect of the self-force, the 
orbital dynamics are computed setting v = in F, as 
well as in A(r) and D(r), i.e., for a Schwarzschild BH. 

Results. Figure Q] shows the absolute value of the 
difference between the Newtonian-normalized (F^ cwt — 
2>2v 2 v w /5) Teukolsky and EOB (uncalibrated and cal- 
ibrated) fluxes as a function of the orbital velocity v. 



The Teukolsky flux includes energy both radiated to in- 
finity and absorbed by the BH's event horizon. The 
uncalibrated Pade-flux (Fj) and p-flux are computed 
through 5.5PN order, but in the Pade flux we also add 
horizon absorption corrections [27| and set v po i e to the 
Schwarzschild light-ring value. The uncalibrated Taylor- 
flux (i.e., the PN Taylor-expanded flux [28[) gives a resid- 
ual about five times worse than the uncalibrated Pade 
and p fluxes. The calibrated Pade-flux (Fg ) is computed 
through 6.5PN order, including the horizon absorption 
corrections, and calibrating v po ie and the 6PN and 6.5PN 
coefficients Tn and J-13; see [28[ for details. The cali- 
brated p-flux is computed through 6PN order, without 
horizon absorption corrections, and calibrating the 6PN 
coefficients Cg 22 in p 2 2 and the 5PN coefficients C5 33 in p 33 ; 
see [ll| for details. The calibration is here performed 
via a least-squares fit to the numerical Teukolsky flux. 
For velocities v S [0.01, 0.1] the agreement is better than 
10~ 8 with a best agreement of 10~ 13 near v — 0.01 for 
all models. 

Comparisons of Teukolsky-based and EOB waveforms 
are performed once they are aligned in time and phase. 
Such an alignment guarantees that the fitting factor 
is maximized over time and phase of coalescence in a 
matched filtering calculation with white noise [13] ■ The 
alignment procedure depends rather sensitively on the 
alignment window chosen. We choose to align the wave- 
forms in the low- frequency regime, i.e., in the inter- 
val t £ [0, 64]A GW , where A GW is the GW wavelength, 
t ~ (0,0.006M) [t ~ (0,0.013M)] months for system-I 
[system-II], to a level of 10 -10 [10 -6 ] rads in the phase 
for system-I [system-II]. We have checked that choosing 
any interval window of width < 2 9 A GW changes the fi- 
nal dephasing by less than 10~ 3 rads and the relative 
amplitude difference by less than 10~ 6 . 

In the left panel of Fig. [2] we plot the absolute value of 
the phase difference, or dephasing, between the dominant 
mode of the Teukolsky-based and EOB waveforms as 
a function of time in units of months. We find that after 
2-years the dephasing is ~ 40 (3000) rads for system- 
I (system-II) when using the EOB-model with Taylor- 
flux (not shown in the figure) [28[ , a result in qualitative 
agreement with previous investigations [2^]. The EOB 
model with uncalibrated Pade-flux at 5.5PN has a de- 
phasing of ~ 5 (530) rads for system-I (system-II), which 
can be reduced to ~ 0.1 (0.01) rads if we employ the 
calibrated Pade-flux at 6.5PN. The EOB model with un- 
calibrated p-flux at 5.5PN has a dephasing of ~ 10 (530) 
rads for system-I (system-II), which can be reduced to 
~ 2 (0.8) rads if we consider the calibrated p-flux at 
6PN. 

In the right panel of Fig. [2 we compare the amplitude 
of the dominant mode A22 = I/122I, computed in the EOB 
and Teukolsky frameworks. After 2-years of evolution, 
both the calibrated Pade- and p-flux EOB models have a 
disagreement of ~ 10 -5 for system-I and ~2x 10~ 3 for 
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P-flux 5.5PN (uncal) 




t [Months] t [Months] 

FIG. 2: Absolute value of the dephasing (left) and fractional amplitude difference (right) of the dominant GW (2, 2) mode as a 
function of time in months. Again, with the introduction of calibrated higher-order terms, the differences are small even over 
a full two year coherent integration. 



system-II. Such a phase and amplitude agreement is fan- 
tastic when one takes into account the 2-year length of 
observation, during which the binary of system-I (system- 
II) evolves over ~ 2 x 10 6 (~ 9 x 10 5 ) rads. Quite in- 
terestingly, we find that if we switch on the relative v 
terms in the 3PN EOB Hamiltonian Eq. ([TJ (conserva- 
tive self- force) and in the flux (dissipative self- force 1 ) the 
dephasing, for the EOB-model with /o-flux at 6PN, in- 
creases to ~ 27 (6) rads for system-I (system-II), while 
the Newtonian normalized amplitude difference increases 
to 4 x (2.5 x 1CT 3 ) for system-I (system-II). We no- 
tice that the main effect comes from the dissipative self- 
force, a result consistent with [3(| for circular orbits (see 
e.g. [aH-ill for more details on the PN self- force). 

We also compare the strongest higher harmonics us- 
ing the EOB model with Pade-flux at 6.5PN. In the case 
of the (I, m) — (3, 3) and (£, m) — (4, 4) modes we find 
dephasings of - 0.14 (0.07) and - 0.18 (0.09) rads, and 
normalized amplitude differences of ~ 6 x 10 -5 (4 x 10 -3 ) 
and ~ 3 x 10~ 4 (9 x 10~ 3 ), for system-I (system-II). 
These dephasings are comparable to those found for the 
(£, m) — (2, 2) mode because in both frameworks the GW 
phase (and frequency) can be computed directly from the 
orbital phase (and frequency), up to errors of less than 
~ 1 rad over a 2-year integration. As a consequence, 
the above comparisons are almost entirely governed by 
the trajectories of the test particle. Finally, we find 
that higher harmonics contribute significantly less to the 
signal-to-noise ratio relative to the (2, 2) mode. In par- 
ticular, we computed the signal-to-noise ratio averaged 
over beam-pattern functions with a noise spectral den- 
sity that includes white-dwarf confusion noise. Including 



up to £ = 5 (£ — 7) for system-I (system-II) guarantees 
a recovery of 97% of the total signal-to-noise ratio, with 
the £ = m modes the most dominant. 

Data Analysis Implications and Discussion. The above 
results have demonstrated that the EOB framework can 
be used to model EMRIs for LISA data analysis pur- 
poses, with the advantage of allowing for the consistent 
inclusion of both dissipative and conservative PN self- 
force terms. In addition, such terms allow the construc- 
tion of waveforms from intermediate-mass ratio inspirals, 
where first-order BH perturbation theory is expected to 
fail. The comparisons made here, however, serve only 
as a proof-of-principle, as one must now generalize the 
formalism to more generic spinning EMRIs, and more 
complicated orbital geometries. 

The EOB framework also allows us to provide, for 
the first time, a metric-based estimate of the number of 
templates needed for EMRI systems in LISA data anal- 
34[. As a coherent 2-year integration in the 
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1 Sometimes all of the energy loss due to radiation, is considered 
part of the dissipative force (even the v = part), but here we 
refer only to the f-dependent terms in the flux. 



ysis 

search of EMRIs is computationally prohibitive, a hier- 
archical search that collects power from coherent searches 
of shorter segments was proposed in [l|. The maximum 
segment length set by computational limits in such a hi- 
erarchical search is estimated to be less than 2 months. 
For a 2-month evolution, we estimate that one requires 
less than 10 7 EOB templates to cover the template bank 
with a minimal match of 0.97 in the total mass range 
(10 5 -10 6 )M Q and mass ratio range (10~ 4 -10~ 5 ). 
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